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ABSTRACT 


This report presents a computer program written in 
Fortran IV for the IBM 360 that is a simulation of the 
illumination on the ground during the descent of an aircraft 
parachute flare from ignition to burn out. The effect of air 
density on the velocity is taken into account by a numerical 
technique. The illumination on horizontal and vertical surfaces 
on the ground are considered. For the surface of interest 
the area consisting of those points having at least a certain 
value of illumination is computed. The program searches for 
the ignition altitude for which this area is maximized over 
the burn time, finds the ignition altitude for which the 
flare burns out at a chosen altitude, or simulates the descent 
with ignition at a chosen altitude. Atmospheric transmission 
is not considered in this report. 

This study was supported by Naval Air Systems Command, 
AIRTASK NO. A35-532/323/69F17-546-502, Work Unit No. 2. 
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I. INTRODUCTION 


This report presents a computer program written in 


Fortran IV for the IBM 260 that is a simulation of the descent 


of an aircraft parachute flare from ignition to burn out. 


The assumptions made in tne sinulation are: 


i 


The flare is considered to be a point source 
of light. 


The candlepower of the flare is constant 
throughout the burning time. 


The burning rate of the flare is constant 
throughout the burning time. 


Vie flares descend vertically in a straignt 
ine. 


The acceleration of gravity is constant and 
equals 32.174 ft/sec. 


Air density variation with altitude above sea 
level is given by 
o(H) = .07513 exp (-.000031582 H) 
where p(H) is the air density in lbs/ft? and H 
is the altitude in feet. 
The drag force on the flare due to air resistance 


is proportional to the square of the velocity, 
i.e., F = KV2, 


The drag constant is not affected by the loss of mass. 


The ground is a flat surface. 


157 
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Two types of surfaces are corsidered. The first is a 
horizontal surface on the ground. The illumination on this 


surface is 


ae 
Ey = pz [0S © 


where I is the candlepower, D is the distance from the flare 
to a point on the ground, and © is the angle between the line 
D and the normal to the ground. 

The second surface is a vertical surface at ground 


level. The illumination on this surface is 


Ey = <. sind, 

The program searches for the ignition altitude which 
results in a burn out at a chosen altitude or which produces 
the maximum area illuminated to at least a chosen value over 
the burn time. A run also can be obtained for any chosen 
ignition altituae. 

A listing of the program and a sample printout is in 


the appendix. 
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II. DERIVATION OF EQUATIONS 
A. Determination of Altitude During Descent 

For the case of variable mass, variable air density, 
and drag proportional tc the square of the velocity, the 
differential equations of motion have no exact solutions. 
Therefore, there are no simple, closed-form mathematical expres- 
sions which can be used to calculate altitude as a function 
of tire for a descending flare system. A numerical technique 
used in this program is a special one-dimensional case of 
the method used and tested by Chipman) in a two-dimensional 
flight. computer program. 

The initial velocity at ignition is calculated as 
the equilhirium velocity (¥) for the given altitude. Equilibrium 
velocity is the velocity where the weight (W) 1s equal to 
the drag force, i.e., W = KV*. Solving for V we get V = VW/k, 
where K = Wk AY The air density, po{H) in lbs/ft°, is a 
function of altitude (H) above sea level in feet and is given by 
o(H) = .07513 exp (-.000031582 H). Cy is the drag coefficient, 
A is the drag area in square feet, and g is the acceleration 
of gravity (32.174 ft/sec). The quantity Cp A/2g is used 
in the computer program as DK (the drag constant). The 


next section describes a method of computing DK. 
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By substituting the expression for K into V = \W/K 


we get the equation for the initial velocity: V = - Vor? 


-—— 


where H is the ignition altitude above sea level. Since 

the positive direction is defined to be up, the descent velocity 
is made a negative quantity. To find the altitude and descent 
velocity during descent a numerical technique is used which 
involves making approximate computations over smal] intervals 

of time. Following is a description of the iteration scheme 
used in the program. 

Let V equal the velocity et the beginning of an interval 
of time At. Two forces, drag in the positive direction and 
weight in the negative direction, are acting to produce an 
acceleration. The total force (F) is -W+KV2. By substituting 


= 2 2 
in the equation a = F/m we get a, = aS or a) = -g + aay 


for the acceleration at the beginning of the interval. 

Letting Z = Kg/W, we have a, = -g + ZV2. Assuming 
this acceleration remains constant during the interval 2 first 
approximation to the velocity at the end of the interval, found 
by using the equation Ve = LF + at, is V; = V+ a,dt. Using 
this velocity the approximate acceleration at the end of 
the interval is a, = -g + ZV,*. A second approximation to the 


velocity at the end of the interval can be found by using 
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the average of the acceleration at the beginning of the 
interval and the approximate acceleration at the end of the 


interval. This gives 


aj+ ao 
V>=V+ a) at 


The average of the velocity at the beginning of the i:.‘erval 
and the second approximation fcr the velocity at the end of tre 
interval is used for the average velocity over the interval. 
Using the equation S = 5, + vt,the altitude at the end of the 
iteration is H = nas (3) At where Ho is the altitude at 
the beginning of the interval. 

Within each interval, the air density, o(H), iS calculated 
based on the altitude at the beginning of the interval and 
the weight is corrected by subtracting the burning rate times 
At. The values at the end of an iteration replace the initial 
values and the process is repeated. 

In the program at is set at 0.1 second. Since for a 


parachute system the descent rate is quite low, this time interval 


is short enough for good approximations. After every 100 iterations 


a printout of time, altitude, and descent velocity is given. 
This gives data for every 10 seconds until the end of the burning 


tim:. Time is equal to zero at ignition. 


ake ee ee ee 
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B. Computation of Drag Constant 
To compute the drag constant DK, solve the equation 


for the initial velocity 
? W 
v= Vote 


for DK. We get 


Pe Ww 
mks SURIVZ 
This equation holds when the system is at its equilibrium velocity 
for a ,iven altitude and weight, which is approximately true 
for a parachute system. Therefore, this equation can be used 
by assuming that the system will be at its average descent velocity 


and at its average altitude when half the composition is burned. 


Example: 
Weight of Flare and Suspension 21 pounds 
Weight of Composition 17 pounds 
Average Descet Rate 7.5 ft./sec. 


Average Altitude (above sea level) 2,000 feet 
The weight when one-half of the composition is burned is 


12.5 pounds. The air density at 2000 feet equals .07053. 


7 W _ 12.5 Ibs. sate . er) 
DK = S(2000)WZ ~*~ TOMES Th. x (75 Ft/secy™ ~ 7:10" Ftrsec’ 
#3 
6 
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C. Area of Illumination for a Vertical Surface 
For a fixed altitude H (in feet) and flare intensity 


I (in candles) the illumination E (in lumens/ft2) 


FIGURE 1 


on a vertical surface at point P on the ground is given by 


* THZERZ 372 (1) 


This equation can be put in the form 


3E\7H*-1¢ 
R& + 3H2R4+ Pape ae R2 +H& = 0 (2) 
V 


which can be considered as; a cubic equation in R2. Choosing 


a value of Ey (say Ey), the condition (see reference 2) for 


two positive real roots is 


(3) 


a 


ee Che eee aes oe 
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If the flare is above this altitude there wil? be no area il- 
luminated to at least Ey and the radii are both set equal to 
zero. For H < Ho the square root of the two positive real 
roots (R; and R») define an area on the ground in which the 


illumination on a vertical surface is greater than or equal to 


Ey: 
Point on ground 
directly below. 
the flare. ~~~: Area in which the 
<=, illumination on a 
vertical surface is 
at least EM 
D. Area of Illumination for a Horizontal Surface 
The illumination on a horizontal surface at pcint P 
is given by 


(H2#R! 3/2 (4) 
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where H is the altitude of the flare. Choosing a value 


of E, (say Ey) and sclving on R we get 
rows 
R= \/(Ey _ 12 (5) 


This is the radius of the area in which a horizontal surface 


has an illumination of at least Ew 


Point on ground SS SA Sg 

directly below SSN — Area in which the 

the flare. SARA illumination on a 
horizontal surface is 


at least EM: 


FIGURE 3 


E. Optimum Ignition Altitude for a Horizontal Surface 
The first step is to find the height for the particular 
flare intensity that will illuminate the maximum area to the 
required minimum illumination. This would be the ideal height 
to burn all of the flare. However, since the flare is descending 
during burning, this is impossible. Therefore, it seems reasonable 


that the ignition altitude should be above the optimum height 


sa af 
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and the burn out altitude should be below the optimum heig! t 
by about the same amount to yield greatest efficiency. Since 
descent rate varies with alt. tude and remaining weight, the 


optimum ignition altitude cannot be computed exactly. 


The illumination Eu at point P from a flare with intensity 


I at point F is Eu = rr cos ©. By substituting cos 0 = H/D, 


we have E,, = IH/D3. Solving for D and squaring we get D2 = (HI/E,,)2/3, 


By substituting this into the equation R =\D2-H2 we get 


(1 \2/3 
R -\(e2) / -H2. To maximize R for a given Ey (say Ey) 
H 


and I we take the derivative of R with respect to H and equate 


to zero (see reference 3). 


2/3 


eae 
AR + (@ rae 


| 
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Solving for ti we get the optimim heignt: 
/y 
_ {i f\2/3f*_ (t/Ey)2/2 (1 /Ey)? /2 
toe YT = APE «Shite 


In the program the approximate ignition altitude is 
computed by adding to the optimum height the product of the 
estimated average descent rate and one-half the burning time. 
The range for the trial ignition altitudes is found by 
multiplying the approximate ignition altitude by 1.25 for the 
upper limit and multiplying by .75 for the lower limit. The 
difference of these two altitudes is divided by 20 to obtain 
the increment between trial ignition altitudes. 

Calculations are made for ignition altitudes at the 
upper limit first. It is possib’e for the altitude to be 
high enough that there is no area illuminated to the minimum 
requirement. The integration of area versus time is done 
by a numerical technique (the trapezoidal rule). The height 
is found at regular time intervals during the descent by 
the method described in section A of part I‘. After every 100 
iterations, the area which is illuminated to the minimum 
illunination is calculated. The value of area-time for each 
interval is computed by multiplying the time interval by the 


average of the area at the beginning of the interval and the 
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area at the end of the irterval. The total integrai of area 

versus time is found by summing over all the burning time. 

The number of iterations and the time of the iterations determine 

the time interval for imiegration. 
The next trial ignition altitude is found by subtracting 

the distance that was calculated for the increment between 

trial ignition altitudes. The same procedure is repeated for 

this trial ignition altitude. When the integral of area over 

the burn time for the present drop becomes less than the integral 

for the previous drop the ignition altitude for the previous 

drop is the estimate the optimum ignition altitude for the 

given illumination level. The true optimum ignition altitude is then 

within the last increment of the estimate. 4 smaller interval 

can be found once the maximum value of the inéegral is past 

and a new set of twenty points computed. | 
In case the flare hits the ground while burning for any of 

the ignition altitudes, no lower ignition altitude is tried | 


because the integral of area-time would only decrease. 


SS ea ee ee 
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III. INPUT TO PROGRAM 


A. Data Format on Cards 


Card FORTRAN Name Columns Format 

1 NPROG 1-2 12 

2 cD 1-10 E10.5 
TW 11-20 F10.0 
WC 21-30 F10.0 
DK 31-40 F10.0 
KOMPUT 41-45 15 
VBAR 46-50 F5.0 
BT 51-55 F5.0 
TD 5€-60 F5.0 
ITER 61-65 15 
HV 66-70 F5.0 
ID} 71-74 A4 
1D2 75-78 A4 
1D3 79-80 A2 

3 H] 1-10 F10.0 

4 EMIN 1-10 F10.0 
HASL 11-20 F10.0 
KODE = 11 


B. Explanation of Terms 
The first data card contains the value of NPROG which 


is a code number used to designate the subroutine to be run. 


NPROG Subroutine 
1 VSB1 (see section IV.B) 
2 VSF1 (see section IV.B) 
3 HSB1 (seé section IV.B) 
4 HSF1 (see section IV.B) 


5 HSO (see section IV.B) 
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The second data card contains the flare parameters. The 
first parameter is the candlepower of the flare (CD). It is 
assumed to be a constant value from ignition to burnout. The 
next parameter is the total we’ght (TW) of the descending system. 
This includes the original amount of composition and the parachute 
system. Tne third parameter is the original amount of composition. 
Botl weights are in pounds. The next parameter is the drag 
constant (DK). Section II.B shows how it is computed. The next 
variable (KOMPUT) provides two ways of reading in DK. If 
KOMPUT = 1 the value read in as DK is used for DK. If 
KOMPUT = 2 the value of DK is computed from the values of VBAR, 
the average velocity of the system (in feet/second), and HV, 
the altitude above sea level] at which the system has an equilibrium 
velocity of VBAR ft/sec. The value of VBAR should be on the 
card when executing the sunroutines VSB1, HSE1, and HSO since 
it is used to comrute an estimate of the ignition altitude. 

The next parameter is the burn time (BT) of the flare measured 
in seconds. The variable TD is the tim increment (in seconds) 
used in the calculation of the height of the flare. The variable 
ITER is the number of iterations used to calculate the height 
of the flare. The product of these two variables gives the 

time between printouts. the three variables IDi, ID2, 1D3 


are used to label the printout with | to 10 characters. 
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The third data card contains the variable Hl. For the 
subroutines VSB1 and HSB1 this value is the ourn out altitude. 
For the other subroutines, VSF1 and HSF1, it is the ignition 
altitude. This card is not required for the subroutine HSO. 

The altitude is measured in feet above the ground (not sea level). 

The height of the ground above sea level (in feet) is 
the second variable on the fourth card. The first variable (EMIN) 
is the minimum illumination level required on the ground measured 
in lumens/sq.ft. The third variable (KODE) takes on integer 
values 0-9 and is used to obtain a run (with the same data) 
uf one of the other subroutines, depending on which subroutine 
was just executed. When KODE takes on a positive value the 


following is indicated. 


NPROG Subroutine to be Run 
1 HSB1 
2 HSF ] 
3 VSB1 
4 VSF 1 
5 VSF1 (at the optimum 


ignition altitude) 
Thus, if the subroutine VSB1 is to be executed first and one 
wishes to obtain a run for the same cata on HSB1, the value of 
KWUE should be a positive integer between 1 and 9 and NPROG get 


to 1. No other subroutine will be executed if KODE = 0. 
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iV. DESCRIPTION OF COMPUTER PROGRAM 
A. Main Procram 
The purpose of the main program is to read in the data 
and to call the indicated subroutines. The first card contains 
NPROG which is an integer used to indicate which subroutine is to 


be executed: 


NPROG Subroutine to be Executed 
1 VSBi 
2 VSF1 
3 HSB] 
4 HSF 1 
5 HSO 


If a blank card is read the program stops. 

The flare parameter card is read in next. If zero 
is read in as the value of ITER, the program will] read in 
another value of NPROG. If ITER has a positive Wale the 
program then computes the drag constant and the burn 
rate of the flare. 

The ignition altitude (if NPROG = 2 or 4) or the 
burn out altitude (if NPROG = 1 or 3) is then read in. This 
card is not needed for HSO (Nrf0G = 5). If zero is read in, 
the program will read in another flare parameter card. 

If a positive value is read in, the program will read in the 
illumination card containing the minirum illumination level, 


the height of the ground above sea level, and the value of 


16 
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KODE. If zero is read in as the value of the illumination level 

the program will read in an ignition (or burn out) altitude 

card. If a positive value is read in, the program will call the 

indicated subroutine. After returning from the called subroutine 
the main program then executes the option indicated by the value 

of KODE. If KODE = 0 another illumination card is read. If 

KODE is a one digit, positive integer the following subroutine 


will be called depending on the value uf NPROG. 


NPROG Subroutine to be Run 
1 HSB] 
2 HSF1 
3 VSBI 
4 VSF1 
5 VSF1 (at the optimum 


ignition altitude) 


Another illumination card is then read. 

The program is set up so that a blank card can be used to 
control which data card is to be read in next. To stop the 
program after executing a set of gata, four blank cards are needed. 
When setting up the data cards, it must be remembered what 


type of data card the program is looking for :.ext. Section V 


shows the data set up for a sample run. 


The value of NPROG is used to indicate which of the 
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B. The Five Main Subroutines 
five main subroutines is to be run by using the following | 

| 


code: 
NPROG Subroutine 
1 VSB] 
2 VSF1 
3 HSB1 
4 HSF1 | 
5 HSO 
The subroutines VSB] and HSB] search for the ignition 


altitude where the flare burns out at a specified altitude above 
the ground. VSB is for a vertical surface and HSB) is for a | 
hroizunta? surface, One can be executed for a set of data just | 
run on the other by use of the variable KODE (See section III). | 
both can determine the ignition altitude by the third trial | ! 
by making a correction to the trial ignition ¢ltitude if the 

burnout altitude is missed by nore than half a foot. Time 1s 

set to zero at ignition. If tne flare hits the ground a cor- 

rection is made to the ignition altitude and the procedure 


started over again. | 
| 
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The subroutines VSF] and HSF1 give a printout for a 
specified ignition altitude and illumination level for a vertical 
and horizontal surface respectively. The variable KODE is used 
to obiain a run on one subroutine with the same set of data 
that was run on the other (see section III). The subroutines 
are ended if the flare hits the ground. 

The subroutine HSO searches for the ignition altitude 
for which the integral of the illuminated area on a horizontal 
surface versus time from ignition to burn out is maximized. 
Section II.E describes the procedure used to obtain the 
estimate of the optimum ignition altitude. An increment 
within which the true optimum ignition altitude lies is also 
given. The variable KODE is used to obtain a run at the 
estimate of the optimum ignition altitude with the same data 


for a vertical surface. 


The subroutine ROOTS determines the radii of the illuminated 


area for a vertical surface (see section II.C.). The subroutine 
HEADER is used to printout data. RHO computes the density of 


the air at a specified altitude above sea level. The subroutine 


RALPH is part of the iteration technique used to compute the altitude 


of the flare (see section II.A. and Chipman‘) , 
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V. SAMPLE RUNS FOR THE MK 45 AIRCRAFT PARACHUTE FLARE 
A. Data Setup 
Data for a sample run is shown in Table 7 on an 
80-column card coding form. The data is set up according to 


the formats listed in section IJI.A. and is as rollows: 


Card 1: NPROG 1 Subroutine VSB] is to be run 


Card 2: CD = 1.65E+06 Candlepower of flare = 1.65x10° 
22.5 Total weight of system (Ibs) 
17.5 Weight of composition (lbs) 
3.14382 Drag Constant 

1 DK is to be the value of the 
drag constant. If KOMPUT = 2, 


the method of section II.B. would 
be used to compute DK. 


' Average descent velocity (ft/sec) 
BT = 180. Burn time (sec) 


TD = | Time increment (sec) 


Number of iterations for the 
computation of the altitude 
of the flare 


HV - Altitude at which the system 
has VBAR as descent velocity. 
Not needed since KOMPUT = 1. 


1D1,1D2, 103 Run identification for 
MK 45-8 printout 


—" 
—j 
Mm 
=) 
" 
~ 
oS 
Lo ] 


Card 3: H] = 300. Burn out a ENE (feet) 
since NPROG = 


Card 4: EMIN = .02 Minimum illumination level 
(Jumens/sq.ft ) 


< za Sz 

o So FR = 

> = 

p=] ~o wt it fl 
(ax 

i —- 

fore] 
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Card 5: 


Card 6: 


Card 7: 


Card 8: 


Card 9: 


Card 10: 


Card 11: 
Card le: 


Card 13: 


HASL = 0. 
KODE = 1] 
EMIN = .05 
EMIN = .10 
Blank 

Blank 

Blank 

NPROG = 5 


Same as Card 2 


EMIN = .02 
HASL = Q. 
KODE = } 
Blank 
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Ground is at sea level 
Run on HSB7 with same data desired 


Another illumination level to 

be run on VSB? for a burn out 
altitude of 300. The resc of the 
card is the same as card 4. 

Hence a run on HSB] is desired. 


Another illumination level. The 
rest of the card is the same as 
cards 4 and 5. 


This blank card indicates that 

there is no more illuminition levels. 
The program will be waiting for a 
burn out altitude card. 


No more burn out altitudes, program 
will be looking for a candle 
parameter card. 


No more candle parameters, program 
will be looking for another value 
of NPROG. 


The optimum ignition altitude 
subroutine (HSO) is to be run 


Minimum illumination level 
Ground at sea level 


A'run on VSF1 with the same data 
is wanted. The ignition altitude 
will be the optimum ignition 
altitude. 


No more illumination levels. 
Program will be expecting an 
ignition (or burn out) altitude. 
Next card should be blank. 
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Card 14 Blank No more values of 1]. The 
program will be looking for 
a candle parameter card. 


Card 15 = Blank No more candie parameters. 
Program will be looking for 
another value of NPROG. 

Card 16 Blank No more subroutines to be run. 
Terminates program. 


B. Description of Printout 


A listing of the Fortran IV computer program written 
for the IBM 360 is given in Appendix A. A sample printout for 
the data setup in section V.A. is given in Appendix B. 

On the printout CD is the flare intensity in candles. 
EMIN is the minimum illuminatiun level in lumens/sq.ft. 

The total weight of the system (TW) is measured in pounds and 
includes the weight of the parachute. The weight of the 
composition (WC) is measured in pounds. It is used along with 
the burn time (BT) in seconds to deternine the burn rate. The 
variable DK is used to determine the velocity of the flare 

and depends on the cross-sectional area of the parachute and 


its drag constant (see section II.A.). HASL is the Feiyht 


of the ground above sea level. 
The time after ignition is given every ten seconds. 
The height of the flare and its velocity is given for ecch point 
in time. The radius (for horizontal surfaces) or tie radii 
(for vertical surfaces) is given as well as the area which has 


the required minimum illumination. The integral of the area 


22 
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illuminated versus time from ignition to each point in time 
is printed out. The average area illuminated since ingition 
is also computed and printed out for each point in time. 

For burn out altitude of 300 feet, runs were obtained 
for minimum illumination levels of 0.02, 0.05, and 0.10 
lumens/sq.ft. on both the vertical and horizontal surface sub- 
routines (see section II). The optimum ignition altitude 
for aminimum illumination level of 0.02 lumens/sq.ft. was obtained 
as well as a printout for the vertical surface subroutine at 


the optimum ignition altitude. 
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TABLE 


DATA PROCESS™NG 80-COLUHN CARD FORMAT 
/ 


DATA Se 7T=ur Far Saree re 


QND-NAIX 5230/2 (6-65) 
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